#!/usr/bin/perl 
if($#ARGV!=3){
    die "syntax: a1 a2 f NIT\n";
}

($a1,$a2,$w10,$NIT)=@ARGV;
$it=0;
$d=$a1-$a2;
$s=$a1+$a2;
$w1=$w10;
$w2=1-$w10;
$filename="a1$a1"."a2$a2"."w1_0$w10".".dat";
open(IN,">$filename");
while($it<$NIT){
    $avgw=($w1+$w2)*0.5;
    $oldw1=$w1;
    $oldw2=$w2;
    $b=0.00001;
    $maxdW=-10000000;
    $bmaxdW=0;
    while($b<=1){
	$z=$b/$d;
	$z2=$z*$z;
	$sqrt1pz2=sqrt(1+$z2);
	$lp=0.5*$d*($s/$d-$z+$sqrt1pz2);
	$lm=0.5*$d*($s/$d-$z-$sqrt1pz2);
	$vpnorminv=1/sqrt($z2+(-1+$sqrt1pz2)*(-1+$sqrt1pz2));
#	$vmnorminv=1/sqrt($z2+(-1-$sqrt1pz2)*(-1-$sqrt1pz2));
#	$v1p=$z*$vpnorminv;
#	$v1m=$z*$vmnorminv;
#	$v2p=(-1+$sqrt1pz2)*$vpnorminv;
#	$v2m=(-1-$sqrt1pz2)*$vmnorminv;
	$v1p=$z;
	$v2p=(-1+$sqrt1pz2);
#	$v1m=$z;
#	$v2m=(-1-$sqrt1pz2);
	$omegap=($w10*$v1p+(1-$w10)*$v2p)*$vpnorminv; #projection of the initial condition on v^\pm
#	$omegam=($w10*$v1m+(1-$w10)*$v2m)*$vmnorminv;
	$f=$omegap*($v1p+$v2p)*$vpnorminv;
#	$fcontrol=$omegam*($v1m+$v2m)*$vmnorminv;
#	print "$b $f $lp $lm $maxdW ", $f*$lp+(1-$f)*$lm," \n";
	if($f*$lp*exp($lp*$it)+(1-$f)*$lm*exp($lm*$it)>$maxdW){
	    $bmaxdW=$b;
	    $maxdW=$f*$lp*exp($lp*$it)+(1-$f)*$lm*exp($lm*$it);
	}
	$b*=1.001
    }    
    $b=$bmaxdW;
    if($it==0){
	$b=1;
    }

    $z=$b/$d;
    $lp=0.5*$d*($s/$d-$z+sqrt(1+$z*$z));
    $lm=0.5*$d*($s/$d-$z-sqrt(1+$z*$z));
    
    $w1=(1+$a1)*$oldw1+$b*($avgw-$w1);
    $w2=(1+$a2)*$oldw2+$b*($avgw-$w2);
    $totw=$w1+$w2;
    $ratio=$w1/$w2;    
    $it++;
    print IN "$it $totw $ratio $w1 $w2 $b $maxdW\n";
}
